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Abstract. We report results of a Wang-Landau study of the random bond square 
Ising model with nearest- (J„„) and next-nearest-neighbor (Jn„„) antifcrromagnetic 
interactions. We consider the case R — Jnn/Jnnn = 1 for which the competitive nature 
of interactions produces a sublattice ordering known as superantiferromagnetism and 
the pure system undergoes a second-order transition with a positive specific heat 
exponent a. For a particular disorder strength we study the effects of bond randomness 
and we find that, while the critical exponents of the correlation length v, magnetization 
f3, and magnetic susceptibility 7 increase when compared to the pure model, the ratios 
(i/v and 7/^ remain unchanged. Thus, the disordered system obeys weak universality 
and hyperscaling similarly to other two-dimensional disordered systems. However, 
the specific heat exhibits an unusually strong saturating behavior which distinguishes 
the present case of competing interactions from other two-dimensional random bond 
systems studied previously. 
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In the last three decades, the effect of quenched randomness to the critical behavior 
of statistical models in two- (2D) and three-dimensions (3D) has been the subject of 
intense studies. First-order transitions are known to be dramatically softened under 
the presence of quenched randomness [H [21 El IH E] , while continuous transitions may 
have their exponents altered under random fields or random bonds pi El [7]. There 
are some very useful phenomenological arguments and some, perturbative in nature, 
theoretical results, pertaining to the occurrence and nature of phase transitions under 
the presence of quenched randomness [S, El El [ini [HI [12J- The most celebrated 
criterion is that suggested by Harris [6]. This criterion relates directly the persistence, 
under random bonds, of the non random behavior to the specific heat exponent 
ap of the corresponding pure system. If is positive, then the disorder will be 
relevant, i.e., under the effect of the disorder, the system will reach a new critical 
behavior. Otherwise, if ap is negative, disorder is irrelevant and the critical behavior 
will not change. The value = is an inconclusive, marginal case. The 2D Ising 
model falls into this category and although it is the most studied case, it is still 
controversial [Ta[H[l5l[l6l[I7in8lE9lEQlEIl[22lE3l[2a[25l[2^ In general 

and despite the intense efforts of the last years on several different models, our current 
understanding of the quenched randomness effects is rather limited and the situation 
appears still unclear for both cases of first- and second-order phase transitions. 

The present Letter is the first investigation of the bond disorder effects on an 
interesting 2D model with competing interactions. We consider the square Ising model 
with nearest- (J„„) and next-nearest-neighbor {Jnnn) antiferromagnetic interactions for 
a certain value of the coupling ratio R = J nn I Jnnn = 1- For this value of i?, the pure 
system undergoes a clear second-order phase transition (from the superantiferromagnetic 
(SAF) state to the paramagnetic state) and accurate estimates of critical exponents have 
recently been reported [281 [29]. Since the value for the critical exponent ap of the specific 
heat of this generalized Ising model is very close (almost identical, see the discussion 
below) to that of the 2D three-state Potts model (a^ = 1/3) our choice of studying 
this case, closely follows the motivation of similar numerical studies performed earlier 
by Kim |30] and Picco [31] on the 2D random bond three-state Potts ferromagnet. In 
other words, due to the difficulties and possible crossover effects in the marginal case 
of Qfp = (2D Ising model) it is desirable to study here the critical behavior induced 
by disorder in a case where the pure model has a positive specific heat exponent and 
according to the Harris criterion [6] is expected to reach a new critical behavior. Our 
results will be therefore profitably compared to those of Kim [3U| and Picco [31j for the 
2D random bond three-state Potts ferromagnet and possible interesting differences may 
reflect aspects that are due to the different microscopic interactions. 

In zero field, the pure system under consideration, is governed by the Hamiltonian: 

T~^p Jnn ^ ^ SiSj -\- Jnnn ^ ^ S^Sj, (1) 

<«J> (id) 

where here both nearest- (J„„) and next-nearest-neighbor (Jnnn) interactions are 
assumed to be positive. It is well-known that the model develops at low temperatures 
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SAF order for R = Jnn/Jnnn > 0.5 ES] and by symmetry the critical behavior 
associated with the SAF ordering is the same under J„„ — > — J„„. We will consider here 
only the case R = Jnn/Jnnn = 1, with J„„ = J = 1. For this case (-R = 1), the system 
undergoes a second-order phase transition, in accordance with the commonly accepted 
scenario for many years of a non-universal critical behavior with exponents depending on 
the coupling ratio R [321 |33l EH ESI ES] . The recent Wang-Landau [37j study of Malakis 
et al [28] has refined earlier estimates [331 ES] for the correlation length exponent v and 
values very close to those of the 2D three-state Potts model z/p(Potts)= 5/6 [38] were 
obtained. From the finite-size scaling (FSS) of the pseudocritical temperatures [28] it 
was found that t'p(SAF;i? = 1)= 0.8330(30) and the subsequent study of Monroe and 
Kim [29j, using the Fisher zeroes of the partition function, yielded a quite matching 
estimate: t'p(SAF;i? = 1)= 0.848(1). Furthermore, from the FSS of the specific 
heat data an estimate for the ratio ap/vp = 0.412(5) was also found [28]. Finally, 
from the magnetic data and in accordance with an earlier conjecture of Binder and 
Landau [33], Malakis et al [28] found additional evidence of the weak universality 
scenario [39] and obtained the values Pp/i'p = 0.125 and 7p/z/p = 1.75. The values 
of the above three ratios of exponents satisfy the Rushbrook relation, assuming that 
Up = 0.8292, which is very close to the estimate obtained from the shift behavior of 
the SAF R = 1 model, thus providing self consistency to the estimation scheme. From 
these results, it is tempting to conjecture that the SAF model with R = 1 obeys the 
same thermal exponents with the 2D three-state Potts model [up = 5/6 = 0.833 . . . and 
ttp = 1/3 = 0.333 . . . [38]), but the respective values of the magnetic critical exponents 
are different {Pp/up = 2/15 = 0.133 ... and 7p/i/p = 26/15 = 1.733 . . . [38]). 

In the present study, we consider a particular type of bond disorder, the same for 
both nearest- and next-nearest-neighbor spins i and j according to the following bimodal 
distribution 

P{J^J) = lm,-Jl)+S{J^,-J2)]■, ^^^ = 1; ^ = f = 0-6- (2) 

The resulting disordered (random bond) version of the Hamiltonian defined in 
equation ([T]) reads now as 

Ti. = JijSiSj + JijSiSj. (3) 

This particular choice of the disorder strength r = 0.6 is strong enough, as will be shown 
below, to observe dramatic saturation effects on the originally diverging behavior of the 
specific heat of the pure model. Only this case will be considered in this Letter; the more 
general case for other values of the disorder strength together with a comparative study 
with the 2D random bond Ising model and also details of our numerical scheme will be 
presented in a longer paper. The present Wang-Landau study [37] closely follows our 
recent implementations of an energy restricted sampling, known as critical minimum 
energy subspace (CrMES) [10] appropriately adapted to the study of systems with 
complicated free-energy landscapes, such as the random-field Ising model ^41j. We 
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Figure 1. Size dependence of tiie maxima of the specific heat for the pure (filled 
squares; data taken from reference [28]) and the random bond model (filled circles). 

impose periodic boundary conditions on square lattices with linear sizes L in the range 
L = 20 — 120 and simulate relatively large ensembles of 100 disorder realizations. 

Each disorder realization is repeatedly simulated up to 4 times with different initial 
conditions. Furthermore, in our implementation of the Wang-Landau repetition process 
thermal properties are calculated at two different Wang-Landau levels and this practice 
enabled us to estimate statistical errors. The statistical errors of the Wang-Landau 
method (WL-errors) used for the estimation of thermal and magnetic properties of 
a particular realization were found much smaller than the statistical errors coming 
from the fact that we used, for disorder averaging, a finite number of 100 realizations. 
Therefore, the WL-errors are not shown in our graphs, since in all cases, they are much 
smaller than the symbol sizes, whereas the latter errors of "finite disorder sampling" 
(fds-errors) are considerable and are presented in all our figures as error bars. The 
mean values over disorder are denoted as [. . .]ai,, the corresponding maxima as [. . 
and the individual maxima as [. . .*]av Since in our fitting attempts we have used mainly 
data from the peaks of the disorder averaged curves (i.e. [C]*„), their fds-errors are the 
relevant statistical errors and have been determined by two similar methods. Using our 
runs, organized in 4 groups of 25 realizations for each lattice size, an application of 
the jackknife method [12] and a straightforward 4-point variance calculation (blocking 
method) [42] were undertaken using the corresponding 4 peaks of the averaged curves, 
for all thermal and magnetic properties studied. It appears that the jackknife method 
yields some reasonably conservative errors, that are about 15 — 20% larger than the 
corresponding calculated standard deviations. These jackknife errors are shown as error 
bars in all our figures and have been used in all our fits. Finally, let us point out that 
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Figure 2. (a) Simultaneous fitting of the pseudocritical temperatures of the average 
specific heat (filled squares) and magnetic susceptibility (filled circles) . (b) Log- log plot 
of the maxima of the average logarithmic derivatives defined in equation ([5|). Linear 
fits are applied for L > 30. Error bars are explained in the text. 



sample-to-sample fluctuations for the individual maxima (such as [x*]ai)) become very 
large with the lattice size. The deflnition of the order-parameter follows reference |28j . 
usine; the four sublattice mae; netizations: M = J2t=i |M|/4. 

Let us start the presentation of our results with the most striking effect of the bond 
randomness on the speciflc heat of the square SAF model. In flgure [U we contrast the 
size dependence of the speciflc heat maxima of the pure (fllled squares) and the random 
bond model (fllled circles). The suppression of the speciflc heat maxima is clear for the 
disordered case, even for the smaller sizes shown and this behavior should be compared 
with the behavior of the speciflc heat of the above mentioned previous studies concerning 
the 2D random bond three-state Potts ferromagnet [301 |3T] . It will then be observed 
from these comparisons that, in full disagreement to our flnding of a strong saturating 
speciflc heat for the SAF R = 1 model, in the case of the 2D random bond three-state 
Potts model one obtains a still diverging behavior for disorder strengths r = 0.9, 0.5, 
and 0.25 [30] and an increasing but progressively saturating behavior is obtained only 
for the very strong disorder r = 0.1 [31]. On the other hand, for the random bond 
SAF R = 1 model, it is evident from flgure [T] that the data of the average speciflc heat 
saturate to a value already from the size of L = 40 and any small variation around this 
value is mainly coming from the fds-errors. Therefore the estimation of the ratio a/i' 
for this model is not possible from the speciflc heat data and the alternative route via 
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the Rushbrook relation will be implemented later. 

In figure[2](a) we present the FSS behavior of two pseudocritical temperatures of the 
model T[z]i^, i.e. the temperatures corresponding to the average specific heat {Z = C: 
filled squares) and the average magnetic susceptibility {Z = x- filled circles). Solid lines 
show a simultaneous fit on both data according to the relation 

Tiz],^ =Tc + hL-"\ (4) 

giving Tc = 1.980(9) for the critical temperature of the disordered model, which is to 
be compared with the corresponding critical temperature Tc;p = 2.0823(17) of the pure 
system [28]. The value for x^/DoF of the above fit, using the jackknife errors, is 0.11. 
Correspondingly, the values of x^/DoF for all our fits vary in the range 0.1 — 0.4. On 
the other hand, using the smaller simple standard deviation errors, one would obtain 
for x^/DoF values in the range 0.2 — 0.7. The above ranges for the ratios x^/DoF reflect 
the goodness of our fits. One could also use the FSS of the pseudocritical temperatures 
defined with the help of the average of the individual maxima of the specific heat and 
susceptibility, i.e. the T[^*]^^, but this choice gives similar results and it is not shown 
here for brevity. A first estimation of the critical exponent v of the correlation length is 
obtained from the above shift behavior and is z/ = 1.080(20), as illustrated in the graph. 
An alternative estimation of the exponent v is attempted now from the FSS analysis of 
the logarithmic derivatives of several powers of the order-parameter with respect to the 
temperature [H US] 

ain(M-) _ (M»E) 

dT "TmT"^^^' 

which scale as L^^'^ with the system size. In figure |2]^b) we consider in double logarithmic 
scale the size dependence of the first- (filled squares), second- (filled circles), and 
fourth-order (filled triangles) maxima of the average over the ensemble of realizations 
logarithmic derivatives. The solid lines shown are corresponding linear fits whose slopes 
provide respectively estimates for 1/v. The estimates in figure [2](b) have an average 
for the correlation length exponent of the order oi v = 1.089. Combining all the above 
estimates we propose an error bound for v of the order of 0.015. Thus, in comparison 
with its value of the pure model, the exponent v for the disordered model shows an 
increase of the order of 30%, reflecting the strong influence of the disorder on the 
thermal properties of the system. It is important to point our here that, our estimate 
is in agreement with the inequality u >2/D derived by Chayes et al [7] for disordered 
systems. 

Turning now the magnetic properties of the model we begin by presenting the 
behavior of the order-parameter at the critical temperature. We present in figure [3]^a) 
in a log-log scale the FSS behavior of the average order-parameter at the estimated 
critical temperature = 1.98. The straight line shows a linear fit for L > 30 with 
a slope of 0.126(5) which is a first manifestation that the ratio (3/v has within error 
bars the value of the pure model. Furthermore, in the inset of panel (a) we plot in a 
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Figure 3. (a) Log-log plot of the average magnetization at the estimated critical 
temperature. The inset shows the log-log plot of the maxima of the average absolute 
order-parameter derivative, (b) Log-lop plot of the size dependence of the maxima of 
the magnetic susceptibilities: [xlav (filled squares) and [x*]at) (filled circles). Sample- 
to-sample fluctuations of [x*]au become much larger than the fds-errors of [x]^^. 



log-log plot the size dependence of the maxima of the average absolute order-parameter 
derivative, defined as 

^^ = {\M\E)-{\M\){E), (6) 

which is expected to scale as L^'^'l^^" with the system size [H US]. Thus, the slope 
of the straight line, which is again a linear fit for L > 30, provides the estimate 
(1 — /5)/z/ = 0.791(11). This estimate when combined with the value for u = 1.089 gives 
a value for j3/v oi the order of 0.127 and using the earlier error bounds for u we propose 
again an error bound of the order of 0.015. Thus, the two estimations are self-consistent 
and our results indicate that although the exponent (3 increases in the disordered case, 
the ratio 13 /v remains unchanged to its pure value, i.e. (ijv = f^p/up = 0.125 |28]. In 
the sequel, we show in panel (b) of figure [3] the behavior of the magnetic susceptibility 
of the model in order to provide estimates for the ratio We present two data 

points: the filled squares refer to the the maxima of the average curve [x]lv^ while the 
filled circles to the average of the individual maxima [x*]ai)- In the latter case the error 
bars shown reflect the sample-to-sample fluctuations, which are, as already pointed out, 
larger than the fds-errors of the corresponding averaged curves and of course much larger 
than the statistical errors. The solid and dotted lines are linear fits for [x]*^ and [x*]av 
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respectively, giving the values ^jv = 1.751(10) and ^jv = 1.757(11), thus providing 
convincing evidence that 7/z/ = 7p/z/p = 1.75 f28], i.e. the ratio 7/z/ maintains the 
value of the pure model. Thus, the ratios [ijv and ^jv for the disordered square SAF 
model appear to be the same with the corresponding ratios of the pure square SAF 
model but different from those of the 2D three-state Potts model. Therefore, our results 
reinforce both the weak universality scenario for the pure SAF model, as first predicted 
by Binder and Landau [33J, as well as the generalized statement of weak universality in 
the presence of bond randomness, given by Kim [30] and concerning also the 2D random 
bond three-state Potts ferromagnet. 

Finally, as discussed above, it is not possible to directly estimate the specific heat 
exponent from FSS of the specific heat data. Yet, having estimated the values for 
/3/z/, 7/z/, and z/, we may estimate a using either the Rushbrook (a + 2/5 + 7 = 2) or 
equivalently, since 2/3/z/ + 7/1/ = 2, the hyperscaling (2 — a = Dv) relation. Both 
relations provide a negative value for the specific heat exponent a = 2 — Dv = 
2 — 2/3 — 7 = —0.173(40), reflecting the early saturation effect. This value of a 
differs significantly from the values estimated for relevant, and even stronger, to ours 
randomness strength for the 2D random bond Potts ferromagnet [30l [31]. One may 
attribute this strongly saturating behavior of the specific heat to the competitive nature 
of interactions which is supposedly responsible for the observed sensitivity of the SAF 
model to bond randomness, since the disorder effects in this case are much more dramatic 
in comparison with the effects observed in other 2D models with simple ferromagnetic 
interactions, such as the 2D three-state Potts ferromagnet [301 [3T]. 

In conclusion, we have applied the Wang-Landau algorithm to investigate the 
interesting effects caused by the presence of quenched bond randomness on the critical 
behavior of the square Ising model with nearest- and next-nearest-neighbor interactions. 
Using standard finite-size scaling techniques, on high accuracy numerical data, we have 
estimated the critical temperature of the disordered model to be well below the value 
of the corresponding pure model and we have extracted values for all critical exponents 
of the random bond square SAF model. These values verify hyperscaling and also 
satisfy the Chayes et al inequality [7] and the weak universahty scenario for disordered 
systems, as stated by Kim [30]. The observed unusual strong saturating behavior of the 
specific heat with a negative exponent a, distinguishes the present case of competing 
interactions from other 2D random bond systems studied previously. 
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